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A numerical method to solve time-dependent Hartree-Fock-Bogoliubov equation is proposed to 
treat the case with Gogny effective interaction. To check the feasibility of the method, it is applied 
to oxygen isotope ^"O and small amplitude oscillations are calculated. The conservation of the 
nucleon numbers as well as energy expectation value is demonstrated. The strength distributions of 
the small amplitude quadrupole oscillations are also shown. 
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I. INTRODUCTION 



The structure and dynamics of exotic nuclei have been 
the main subject of investigation both in the theoretical 
and the experimental nuclear physics. The peculiar fea- 
ture of the exotic nuclei is that the Fermi level is located 
near the continuum levels. In this situation, the nucleons 
near the Fermi level are easily brought to the continuum 
states by pairing correlations. The low-energy excitation 
modes are often in the continuum energy region. 

The general mean-field method to treat the pairing cor- 
relations as well as the mean-field is the Hartree-Fock 
Bogoliubov (HFB) method, which has played a central 
role in the investigation of the static ground state prop- 
erties of the nuclei in a wide area of the nuclear chart 
The typical method to study the excited collective 
states of nuclei on top of the mean-field ground state with 
the pairing correlations is the quasiparticle random phase 
approximation (QRPA) 0, S i, 0, S i, US] ■ 

In the practical QRPA calculations, two major meth- 
ods so far developed are the Green's function method 
through the response function formalism and the diago- 
nalization method of the QRPA matrix. From the the- 
oretical point of view, the QRPA is a small amplitude 
approximation of the time-dependent Hartree-Fock Bo- 
goliubov (TDHFB) equation [ll|. Considering the rich 
results of the time-dependent Hartree-Fock (TDHF) and 
the random phase approximation (RPA) as a small am- 
plitude approximation of the TDHF it is worth- 
while establishing a self-consistent TDHFB framework 
together with the practical method of solving the TD- 
HFB equation. 

The widely used effective interactions both in the 
mean-field and the QRPA calculations are the Skyrme 
interactions with zero-range pairing part or the Gogny 
interaction which consists of both the finite-range parts 
and the zero-range ones. 

In the case of the HFB calculations with the Skyrme in- 
teractions (Skyrme HFB), the Skyrme interaction is used 
for the particle-hole channel, while the pairing interac- 
tion is introduced only for the particle-particle channel. 
Since the zero-range interaction is assumed in the Skyrme 
HFB, it is necessary to set the appropriate cut-off energy 



and choose the optimum parameter set in the pairing 
part [Hi. 

In contrast with the Skyrme HFB, in the HFB cal- 
culations with the Gogny interaction (Gogny HFB), the 
particle-hole channel and the particle-particle channel are 
treated on an equal footing. Therefore, the Gogny inter- 
action is suitable for the formulation of the self-consistent 
TDHFB together with a practical numerical method of 
integrating the TDHFB equation. 

In this paper, we report the formulation of the self- 
consistent TDHFB with the Gogny interaction and a 
numerical method to integrate the TDHFB equations. 
With the aim at illustrating the feasibility of the method, 
we apply the method to the case of oxygen isotope ^"O, 
carrying out the numerical integration of the TDHFB 
equations. 

This paper consists of the following sections: In Sec. 
II, a derivation of the self-consistent TDHFB equation 
is given together with a numerical method of solving the 
TDHFB equation. In Sec. HI, the results of applying the 
TDHFB equation to the oxygen isotope ^"O are shown. 
Section IV is for summary and concluding remarks. 



II. TDHFB EQUATION AND NUMERICAL 
SOLUTION 



1. TDHFB equation 

Let us assume that the Hamiltonian is written as a sum 
of the kinetic energy and the two-particle interactions 
with particle creation (annihilation) operator C^^ (Cq), 



Q/3 



where Tap is the kinetic energy matrix element and Vap-yS 
is an antisymmetrized two-body matrix element. 

In the formal presentation of the TDHFB equation, it 
is convenient to start with the generalized density matrix 



p 



l-p* 



(2) 



2 



where p and k are normal density and pairing tensor, 
respectively, 



r*TjT\ 



(3) 
(4) 



The matrices Uak and Vak are introduced to connect the 
particle operators {C|, Co) with the quasiparticle oper- 
ators /3fe} as follows: 



a 



(5) 
(6) 



The equation of motion of the generalized density ma- 
trix TZ in Eq. ^ is given by 

ihil ^ [H, 7^] , (7) 
with the HFB Hamiltonian Ti given by 

h A 



n = 



-A* -h* 



(8) 



Here, mean field Hamiltonian h and pairing mean field 
A are introduced through the following relations: 



lap 



Ta/3 + Fq/j, Ta/3 = ^ '^a-ypSPS-f, (9) 



(10) 



7(5 



From the definitions in Eqs. ([3]) and (H]), the general- 
ized density matrix TZ in Eq. ([2]) is rewritten as 



n 



V* 

u* 



(11) 



where T stands for transposed matrix. Combining the 
formal solution TZ{t) of the equation of motion in Eq. ([7]) 
given as, 

7e(t) =e-^/'''^^M7^(0)e^^'^^^(^^ (12) 

with the rewritten form of the generalized density matrix 
TZ in Eq. pT|) , we have a formal solution of the matrices 
U and V as follows: 



u*{t) y I U*{0) 



(13) 



The formal solution Eq. (|13p of the matrices U and V is 
equivalent to the equation of motion of the matrices U 
and V given as 

dt I u*{t) J \ u*{t) j ' 



( uit) \ u{t) 
dt [ v{t) J - ^[ v{t) 



(14) 



where the definition of the HFB Hamiltonian Ti in Eq. 
(IHl) is used. This form of the TDHFB equation was used 
by Bulgac in relation with Berry's phase jl4i] . 



2. Numerical method of solution 

The TDHFB equation takes a simple form, being 
similar to the TDHF equation with a TDHF Hamiltonian 
hroHF, 



d 

«fi^V'j(x,t) = hTDHFijj{x,t), 



(15) 



for the wave functions t){j = 1, 2, • • • , N) of N or- 

bitals. The solution of the TDHF equation is calculated 
by making use of the relation, 

V'j(x,t„+i) =e-^^"'"+'''V(x,i„), (16) 

at every time step from t„ to = tn + At with an ade- 
quate Hamiltonian to conserve the total energy 
lE^I. Applying the method in ^ to the TDHFB equa- 
tion (fH|) . we get the solution of the TDHFB equation 
(fT4)l in the form given as 



(n+l) 



At, 



= exp -iili7^(»+i/2) 



(n) 



, (17) 



with an adequate TDHFB Hamiltonan 7^^"+^/^^ at every 
time step from t„ to = i„ + At. 

In the present case with the Hamiltonian ([T|), the ex- 
pectation value of the Hamiltonian ([T]) with respect to 
HFB state |$) is given as 



E = 

a(3 



+ -Taf3Pfla + -K*^pAal3, (18) 



where the mean potential F and mean pairing potential 
A are defined in Eqs. dU and HUl). Since in the TDHFB 
calculation the energy conservation is one of the most 
important conditions to be fulfilled, let us see how the 
energy in (jlSp is conserved with respect to a small vari- 
ation in the matrices U and V. Within the first order of 
the parameter A = ^ the matrices U and V at a, time t 
is changed into the new matrices U' and V' according to 
the following relation. 



U' 
V 



-iX 



h 

-A* 



A 

-h* 



U ~i\{hU + AV) 
V + i\{h*V + A*U) 



(19) 



Here, for the ease of the discussion, let us assume that 
the time increment At is small enough so that we can 
identify the mean field Hamiltonian and mean 

pairing potential A("+i/2) in the TDHFB Hamiltonian 
-^(n+i/2) .^yjth h and A, respectively, at the time t. Using 
the relations in Eq. ([TO)) , the variations in the density p in 
Eq. ^ and pairing tensor k in Eq. Q are represented, 
respectively, as 



3 



p' = V'*V'T 

= (y* - 1\ (AC/* + w*)) {y^ + 1\ (t/^At + y^/it)) 

= p-iA[/i,p] -iA(-AK* + kA*), (20) 
- (y* - zA (AC/* + /iV^*)) (C/^ - i\ {U^h^ + V^A^)) 

= K-i\{AU*U^ + hK + Kh* - pA) . (21) 



Putting these expressions of the density and pairing ten- 
sor up to the first order in the parameter A into the ex- 



where the relations h* = hF , = —k, and A-^ = — A 
are used. The notation Tr stands for taking the trace of 
the matrices. From Eq. ([H]), we see that we can integrate 
the TDHFB equation (fT4l) with a conserved energy ex- 
pressed as in Eq. (fT8|) . setting the time increment At and 
intermediate TDHFB Hamilonian 'H^'^+^Z^') adequately. 

When the Gogny interaction is used, there is a part 
in the mean potential F which comes from the density- 
dependent term through the variation of the density ma- 
trix p, just as in Eq. ([20|l . Owing to the parameter set 
of the Gogny interaction, there is no contribution of the 
density-dependent term to the pairing energy. The con- 
tribution of the density-dependent term is included only 
in the mean-field Hamiltonian h in the HFB Hamiltonian 
([8]). Then, the energy conservation relation ((22l) holds 
when the Gogny interaction is adopted in the Hamilto- 
nian 

In the numerical calculation, we expand the exponen- 
tial function in Eq. (|17p in terms of the power series up 
to the tenth order in the time increment At. The inter- 
mediate TDHFB Hamiltonian is made by us- 
ing the predictor-corrector method at each time step: In 
the predictor-corrector method, the predictor solutions 
U' and V' are calculated according to the method in 
Eq. ^ by using the TDHFB Hamiltonian H^") with 
the quantities h and A, which are made by using den- 
sity p^"-* and pairing tensor k*^"^ at the time respec- 
tively. Then, using the predictor density p' = V'V^ 
and pairing tensor k! — V'*U''^, the intermediate den- 
sity = (p(")+p')/2 and the pairing tensor 



I 

pression of the energy in Eq. (|18p , we have the variation 
of the energy as follows: 



(22) 

I 

^(n+i/2) _ (^ni"^) + k') /2 are made, which enter the 
intermediate TDHFB Hamiltonian n^''+^/'^\ The cor- 
rector solutions U" and V" are calculated according to 
the method in Eq. (fT7|) with the intermediate TDHFB 
Hamiltonian 7i("+^/-^) . From the ideal point of view, this 
process is repeated until the energy is conserved within a 
desired order. In the practical calculations, however, we 
stop the predictor-corrector iterations after the first two 
iterations to save the cpu time. 

The initial condition we adopt in the present calcu- 
lation is of the impulse type: The static HFB solution 
Uq and Vb are changed into the initial matrices C/'-*'-' and 
y^"^ by the relations given by 



= exp i^eQ) Vo = "j2 ^^^^o, (23) 
i/=i 

C/W = exp(-z£Q*)C/o- J2 ^ , Uo, 

(24) 

where the expression Q stands for matrix representation 
of a multipole operator with respect to the numerical 
basis states. In the expression of the initial conditions 
([23|) and (p4|) . the exponential function is expanded into 
the power series with respect to the parameter e up to 
the A^ma^-th order. 



5E = Tr{-a/i[/i,p] -a/i(-AK* + kA*)} 



2 



Tr{-A*A - A*hK - A*Kh* + A*Ap* + A*pA} 



-yTr{ A*A + /i*K*A + K*/iA- A*pA-p*A*A} 
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FIG. 1: Time dependence of expectation value of quadrupole 
moment. 



FIG. 3: Deviation of neutron (proton) number expectation 
value from accurate number 12 (8) in quadrupole oscillation 
in Fig. □ 
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FIG. 2: Total energy vs time in quadrupole oscillation in Fig. 
[U Line labelled with "ground state" stands for the energy 
of the HFB ground state. Curve labelled with "TDHFB" is 
for the energy expectation value in the course of quadrupole 
oscillation. 



III. APPLICATIONS 

With the purpose of studying the feasibility of the TD- 
HFB equation and the method of solution proposed in 
the previous section, we apply the method to the oxygen 
isotope ^"0. We adopt the Gogny force [ii,[i3| with Dl 
parameter set in the two-particle interaction part in the 
Hamiltonian (1). 

As the numerical basis, we make use of the three- 
dimensional harmonic oscillator eigenstates. The two- 
particle matrix elements Vap-yS in the Hamiltonian (1) are 
calculated after the method which was used by Girod and 
Grammaticos 18]. In the present calculations, we set the 
space of basis states so that the relation Ux + ny + Uz < 
N shell = 4 is satisfied, where rix [uy, Uz) is the number 
of quanta of the harmonic oscillator basis states in the 
X (y, z) direction, respectively. The angular frequency 
parameters ujx, ujy, and uJz of the harmonic oscillator ba- 
sis states are optimized under the sphericity condition 
ujx = LiJy = oJz = ^0 so that the HFB energy should be 



FIG. 4: Quadrupole strength distribution of quadrupole os- 
cillation in Fig. [1] Artificial width of 0.6 MeV is used. 

minimum. We have set fiWQ — 16.0 MeV. 

In the initial conditions ([25)1 and ([M]) . we take the 
multipole operator Q to be a quadrupole operator Qap = 
(2z'^ — — y^) expressed as a matrix with respect to 
the basis states labelled by a and /3. The parameter e is 
put to be 1.0 X 10~^, which is small enough so that the 
linearity of the oscillation with respect to e is satisfied. 
The power series expansions of the exponential functions 
in Eqs. ((23)) and (|24)) are calculated up to the N^ax-^^ 
order with N^ax put to be ten. 

In the present calculations, we have omitted the 
Coulomb part in the two-body interactions, which leads 
to shorter cpu time. For the moment, it might not be a 
draw-back to neglect the Coulomb part with the aim at 
studying the feasibility of the method under considera- 
tion. 

In Fig. [l] we display the time variation of the mass 
quadrupole moment (2z^ — — y^) of ^°0. The time 
increment cAi used in the calculation is 0.2 fm, and total 
time step is 19000. After the initial impulse, we can see 
regular small-amplitude oscillations take place. 

In Fig. [51 the total energy in the course of the oscil- 
lation is shown together with the ground state energy of 
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the static HFB calculation. The excitation energy 0.2 
MeV is kept to a good extent in the integration process. 

In Fig. [31 the deviation of the expectation values of the 
nucleon number from the accurate values (8 protons and 
12 neutrons) is displayed with respect to time. In ^'^O, 
the protons arc in the normal state, whereas the neutrons 
are in the superconducting state. Therefore, the integra- 
tion of the equations of motion of the proton orbitals is 
equivalent to the TDHF case, where the occupation prob- 
ability of each one of the orbitals is exactly 1 or 0. Then, 
the total number of protons are conserved within 10^^^. 
The neutron number in Fig. [3] is kept up to around 10~^ 
in the present integration process. This result illustrates 
that we can keep the unitarity of the time development 
operator in Eq. ()17p within a practically satisfying order 
with the following set of parameters such as Nsheii = 4 
and cAt = 0.2 fni in the case of the excitation energy 0.2 
MeV. 

In Fig. m we display the strength function which is 
calculated by the Fourier transformation from the time 
series of the expectation value of the quadrupole operator 
in Fig. [T] The locations of the main peaks in Fig. 2] 
are similar to the results in Ref. 6. The low-energy 2+ 
levels of are located at 1.7 MeV, 4 MeV, 5 MeV, 
and 10 MeV [19j. The three low-energy peaks in Fig. 
m are expected to correspond to these levels. Since the 
space of the basis states is not large enough, the energies 
corresponding to the peaks in Fig. |4]are somewhat higher 
than those of the observed levels. 

On the other hand, the ratio of the lowest energy peak 
around 3 MeV to the high energy one around 24 MeV 
is quite different from the result in Ref. 6. Considering 
that the effect of the continuum is included in the RPA 
calculations in Ref. 6, and our space in the present cal- 
culations is limited within the major quantum number 
Nsheii = 4, the differences of the low-energy peaks might 
come from the space of states taken in the calculations. 



integration method which is widely used in the TDHF, 
i.e., a power series expansion of the time displacement op- 
erator at each time step. Adopting the Gogny interaction 
in the two-particle interaction parts in the Hamiltonian, 
we carried out the numerical calculations of the TDHFB 
equation in the case of oxygen isotope ^"O. 

The accuracy of conservation of the excitation energy 
and the particle numbers is illustrated in the calcula- 
tion of small-amplitude quadrupole oscillation, which is 
started from an impulse-type initial condition. 

The strength function of the quadrupole oscillation is 
calculated. The energies of the peaks of the strength 
function are similar to and somewhat higher than the 
experimental results of the low-energy 2+ levels. 

The relative ratios of the peak heights, on the other 
hand, seem not to be enough to discuss their physical 
contents, in comparison with the results given by the 
QRPA with continuum space included Q. One major 
reason of the situation might be the space of basis states 
in the present calculations, which is not large enough 
with the maximum major shell quantum number Nsheii 
= 4. It is the foremost task to check the convergence of 
the results with respect to the cut-off Nsheii of the space 
of basis states. 

Since the basis states in the present calculations are the 
three-dimensional harmonic oscillator states, the method 
proposed in this article seems to be useful in describing 
the excitational modes in deformed nuclei. The calcula- 
tions in some deformed nuclei are now in progress. 
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